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RESUME 


Dans  ce  rapport >  on  aSveloppe 
d'algorithmes  congus  pour  permettre  a  une 
recherche  de  poursuivre  une  source  de 
droite.  Les  unites  de  recherche  mesurent 
apparente  afin  d'estimer  la  distance, 
composantes  de  la  vitesse  de  la  source, 
gisement  et  de  frequence  sont  sujettes  b 


certains  fon&ements  theoriqves 
ou  plusieurs  unites  mobiles  de 
bruit  qui  se  d^place  en  ligne 
le  gisement  et  la  frequence 
la  fr€quence  r^elle  et  les 
On  presume  que  les  mesures  de 
des  erreurs  alSatoires.  (NC) 
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ABSTRACT 


r  JP.  r^r-=  »?!=?»: 

sound  source  moving  m  a  straight  Une.  frequency  and  velocity 

and  received  frequency  to  infer  range  the  true  frequency  to 

components.  The  measurements  of  bearing  and  frequency  are  0 

random  errors.  (U) 
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1.0  INTRODUCTION 

This  report  presents  a  basic  numerical  treatment  of  the 
two-dimensional  problem  posed  by  the  passive  tracking  of  a  sound  source 
by  one  or  more  searchers  that  measure  the  bearing  and  freouencv  of  the 
source.  Th 2  question  is:  can  the  position,  the  course,  the  speed,  and 
the  frequency  of  the  source  be  determined  with  usable  accuracy  when  the 
observations  about  its  bearings  and  its  freauency  are  made  over  a  period 
of  time  and  in  the  presence  of  error  noise? 

Constant  source  velocity,  i.e.  constant  speed  and  direction,  is 
the  only  case  considered  in  this  report.  Further,  the  measurement 
errors  of  both  bearing  and  frequency  are  assumed  to  be  normally 
distributed,  with  specified  means  and  standard  deviations  and  to  be 
independent  from  measurement  to  measurement. 

Clearly,  complete  tracking  cannot  always  te  accomplished  on  the 
basis  of  freauency  observations  alone,  because  a  source  moving  alone  a 
tangent  to  a  circle  around  a  single  searcher  will  always  exhibit  the 
same  Doppler  effect,  no  matter  which  tangent  this  source  follows.  Of 
course,  the  use  of  multiple  searchers  or  chanees  in  searcher  velocity 
during  measurement  may  yield  additional  information,  but  the  algorithms 
discussed  in  this  report  are  based  on  a  single  fundamental  freauency 
relation,  which  uses  both  bearing  and  freouency  information. 

There  seems  to  be  some  confusion  about  the  exact  form  of  the 
Doppler  equation,  so  I  first  related  the  usual  approximation  for  Doppler 
frequency  to  a  correct  relation  for  observed  freouencv  when  both  the 
source  and  the  searcher  are  moving  relative  to  the  propagation  medium;  I 
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then  numerically  explored  the  accuracy  of  the  approximation.  This*  is 
reported  in  Section  2.0. 

In  Section  3.0,  the  acceptable  approximation  relation  found  in 

Section  2.0  is  used  to  develop  a  fundamental  system  of  linear  eauations 

between  the  source  parameters  and  the  measured  data.  Each  measurement 

provides  one  equation,  and  it  is  assumed  that  the  total  number  of 

measurements,  N,  exceeds  three,  which  is  the  number  of  parameters  to  be 

determined.  The  solution  of  the  resulting  system  of  equations,  in  least 

.  * 

squares  sense,  for  f,  XT,  YT,  is  discussed  and  the  effects  of 

constraints  on  the  region  of  acceptable  solutions  are  considered  in 
Section  3.0.  Since  the  linear  equations  of  Section  3.0  do  not 

explicitly  contain  range,  their  solution  produces  no  Information  about 
this  parameter. 

Several  ways  to  introduce  range  as  a  parameter  in  the  linear 
system  are  discussed  in  Sec.  4.0.  The  most  informative  one  seems  to  be 
the  straightforward  grid  search  for  a  least  squares  minimum  deviation 
between  predicted  and  observed  freauency-only  data.  However,  the  grid 
search,  which  uses  an  internal  iteration,  is  very  time  consuming,  and 
can  be  taken  only  as  indicative  of  possible  information  in  the  data. 
This  information  might  then  be  better  extracted  by  some  other  procedure. 
In  Section  5.0  comments  are  offered  on  sample  calculations  made  with  the 
algorithms  developed  earlier.  The  actual  results  of  those  calculations 
are  not  given.  Iv,  is  worth  commenting  that  the  problem  has  a  unique 
solution  even  on  the  basis  of  only  three  measurements,  provided  these 
measurements  are  free  of  errors.  Of  course,  there  are  time-honored 
methods  for  passive  tracking  using  bearings  only,  but  they  require  more 
observations,  multiple  observers,  or  maneuvers  by  the  observers.  The 
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proof  that  the  Joint  bearing-freaueney  problem  has  a  unioue  solution  for 
three  observations  is  given  in  Aopendix  B. 

In  Appendix  A  the  APL  implementations  of  the  olgorlthms  outlined 
in  the  text  are  detailed .  Much  improved  imolementations  could  be 
produced,  but  those  given  are  adequate  to  test  the  algorithms.  What  is 
really  needed  is  a  good  algorithm  for  determining  range  when  this 
parameter  is  not  a\ailable  as  a  primary  variable. 

The  text  also  includes  some  comments  on  the  effect  of  providing 
estimates  of  one  or  more  parameter  values  to  decrease  the  breadth  of  the 
search . 

This  work  was  performed  at  DREV  in  the  first  half  of  1980,  under 
PCN  32D37,  Tactical  Towed  Array  Study. 
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2.0  THE-EFEQUEIO  SHIFT  RELATION.  AND  ITS  APPROXIMATION 

When  both  the  souroe  of  a  sound  and  the  receiver  are  in  motion 
with  rospeot  to  the  propagation  medium,  a  component  of  the  resultant 
frequency  shift  arises  from  each  of  the  two  motions.  In  the  following 
analysis,  the  absolute  frequency  of  the  souroe  is  represented  by  the 
aooonpanying  absolute  time  between  wave  crests,  or  pulses,  denoted  by  A. 
Frequency  is  the  reciprocal  of  A.  Thus,  absolute  percentage  errors  in  A 
and  in  frequency  are  the  same,  since 

f  «  1/A 
Df  -  -DA/A2 
Df/f  -  -DA/A 

In  Fig.  1,  T  and  T'  are  the  positions  of  the  source,  a  time  A 

apart,  and  S  and  S'  those  of  the  reoeiver  when  the  pulses  emitted  at  T 

and  T'  are  received.  The  time  interval  between  the  reception  of  the  two 

pulses  is  denoted  A  .  The  rest  of  the  notation  used  in  the  diagram  is 

a 

defined  as  follows: 

*  the  vector  velocity  of  the  source  with  respect  to  the  medium 
s  the  speed  of  the  source,  or  magnitude  of  VT 
a  the  speed  of  sound  in  the  medium 

=  the  veotor  velocity  of  the  receiver  with  respect  to  the 
med  ium 

=  the  speed  of  the  receiver,  or  magnitude  of  Vg 
=  the  vector  from  S  to  T,  l.e.  the  reverse  of  the  propagation 
path  for  the  first  pulse 
s  the  distance  from  S  to  T 


Figure  1  -  Geometry  of  'bearing  and  frequency  measurements 


Rg  s  the  veotor  from  S'  to  T',  i.e.  the  reverse  propagation  path 
for  the  second  pulse 
R2  *  the  distanoe  from  S'  to  T' 


I 


! 

j 


|  A  and  B  are  angles  measured  clockwise  from  the  heading  vectors  to  the 

propagation  paths,  as  illustrated. 

i 

r 

To  compare  the  time  interval  between  the  reception  of  the  pulses, 

I  A. ,  with  that  between  their  emission,  A,  we  write  an  equation  statin* 

f  * 

!  *  that  time  along  path  T  *  T'  ■*  S'  must  equal  time  alon*  path  T  ♦  S  ♦  S'. 

i 

t  ♦  Rg/C  »  ♦  R^/C  [1] 

[ 

I 


If  we  introduoe  coordinates  <XT,  YT>  and  <XS,  Ys>  to  desoribe 
positions  of  the  points  T  and  S,  and  let 
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then,  of  course, 

R1  -  (X2  +  Y2)* 


[3] 


and 


((X  +  AXp  -  LXS)2  +  (Y 


«*  ayt  - 


Vs)£)i 


[k] 


where  the  dots  indicate  time  derivatives,  or  velocity  components.  A 
direot  substitution  of  [3]  and  [4]  into  [1]  yields: 


A 

a 


**  A  +  ~  (Rg  -  Rl) 

-  A  +  |  ({(X  +  AX^Xg^+CY  +  AYT-AaYs)2)i-(X2+Y2)i) 


[5] 


Equation  [5]  is  an  implicit  relationship  for  the  unknown,  A  ,  in  terms 

S 

of  the  other  quantities,  all  of  which  are  known.  This  equation  cannot 

be  easily  solved  directly  for  A  ,  but  that  step  is  unnecessary  since  [*5] 

can  be  used  as  an  iteration  equation  which  converges  very  rapidly  to  the 

value  of  A  ,  starting  from  A  =  A. 

8  a 


The  APL  proeram  DEXACT  (Appendix  A)  is  an  iteration  which  serves 

to  calculate  A  to  accuracy  10-1U. 
a 


The  iteration  based  upon  [5]  converses  to  the  true  value  of  A. 

a 

but  it  is  not  in  the  usual  form  of  the  Doppler  equation,  and  a 
closed-form  approximation  which  provides  sufficient  ecouracy  would  be 
more  useful.  To  derive  the  best  linear  approximation  to  [5],  one 

regards  the  quantity  in  brackets,  F2  “  Rj,  for  fixed  A  and  Ag,  as  a 

*  »  »  « 

function  of  the  six  variables  X,Y,XT,YT,XS,YS>  As  such,  it  can  be 
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expanded  in  a  six-dimensional  Taylor  expansion  about  point  X,Y»0, 0,0,0. 
If  this  is  done,  and  if  only  the  first-order  (linear)  terms  are 
retained,  the  approximated  equality  [6]  results: 


A 

a 


4  +  C  (f  4  *T  *  R  4  *T  ‘  I  4a  *S 


R  4a  YS> 


[6] 


Noting  that  X/R  and  Y/R  are  the  components  of  the  unit 
vector  "r^  8  R^/R^  the  dot  or  soaker  product  of  two  vectors  oan  be  used 
to  give  [6]  a  more  symmetric  form. 


A 

a 

A 


1 


1 


+  K 


L  i->- 
C  ri 


[7] 


The  Doppler  shift  is  usually  represented  by  [7). 


The  dot  products  in  [7]  can  be  written  in  terms  of  angles  A  and  B: 


VT  cos  (»  -  A)  s  - V T  cos  A 
Vs  cos  (2w  -  B)  =  Vg  oos  B 


and  [7]  then  becomes 


A 

a 


1  ”  C  VT  C0S  A 


[8] 
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A  better  approximation  to  eouation  [5]  can  be  obtained  with  the 
expression  generally  used  for  Doppler  shift,  particularly  in  the  case  of 
one-dimensional  motion: 


[9] 


where  R  is  the  time  derivative  of  the  distance  between  the  receiver  and 
the  source.  From  Fig.  1,  it  can  be  seen  that 

R  «  V„  cos  (it  -  A)  -  V0  cos  (2t-  E) 

1  o 

«  -V  cos  A  -  V  cos  B 
T  o 

and  hence,  [9]  can  be  written 

A 

~  -  1  -  ~  (VT  cos  A  +  Vs  cos  B)  [9'] 


But  [9']  also  results  from  [8]  because  (Vs/C)cos  B  is  small  with  respect 
to  1;  hence, 


(1  +  |  Vg  cos  B)’1  ■  1  -  £  Vg  cos  B 


approximately.  If  the  term  involving  1/C  is  discarded  in 


~  ■  (1  -  £  VT  cos  A)  (1  -  J  Vg  cos  B) 


[10] 


The  errors  inherent  in  these  approximations  were  sampled  by  the 
APL  function  DAPPROX  (see  Appendix  A)  for  two  signal  frequencies  nearly 
two  octaves  apart.  The  results  were  not  much  different  for  the  two 

frequencies.  The  errors  of  approximation  of  [5]  by  [8]  and  [9]  were  of 

.5  _3 

the  order  of  10  percent  and  10  percent,  respectively. 


It  should  be  stressed  that  tbest  errors  arise  from  an 
approximation  formula  and  are  therefore  not  only  small,  but  also 
systematic.  They  are  not  random  errors.  For  frequency  tracking,  where 
the  true  frequency  is  unknown,  systematic  errors  do  not  influence  the 
accuracy  so  long  as  they  do  not  change  significantly  while  tracking  data 
are  gathered.  Further  numerical  investigation  of  the  errors  involved  in 
using  [7]  to  compute  the  Doppler  shift  reveals  that  they  are  systematic 
and  slowly  varying.  Thus,  it  is  concluded  that  equation  [7]: 


A 

a 

A 


.1 

j 

l 

\ 

\ 


1 
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3.0  THE  FUNDAMENTAL  LINEAR  SYSTEM 

The  basic  eauation  for  tracking  a  sound-emitting  source  by 
measuring  its  bearing  and  frequency  continuously,  or  at  several  discrete 
points  in  time  is: 


f  c  +■?  •  V,. 

a  _  _ S 

f  c+*?.vT 


[11] 


as  derived  in  the  preceding  section,  where: 


f  =  apparent  freauency  to  the  listener 
f  s  true  frequency 
C  =  speed  of  sound  in  the  medium 

r  =  a  unit  vector  pointing  along  the  reverse  of  the 
pulse  propagation  path,  from  searcher  to  source 
V  =  velocity  vector  of  the  searcher,  relative  to  the  medium 

O 

VT  =  velocity  vector  of  the  source,  relative  to  the  medium. 

It  should  be  noted  that,  because"?’  is  a  unit  vector,  the  range 
from  searcher  to  source  does  not  explicitly  appear  in  [11].  Because  of 
this,  range  is  a  secondary  parameter  in  the  search  problem,  and  is 
difficult  to  determine  since  the  measured  quantities  of  freouency  and 
bearing  are  not  sensitive  to  range  variation.  In  the  search  algorithms 
proposed  here,  range  is  not  initially  determined,  but  it  is  left  to  be 
determined  by  a  follow-up  algorithm. 
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If  B  is  the  value  of  the  bearing  angle,  measured  clockwise  from 
the  north,  then  the  unit  vector  r  has  components 

<  cos  B,  sin  B  >  [12] 

In  this  treatment  the  positive  X-axis  is  oriented  toward  the 
north  and  the  positive  Y-axis  toward  the  east.  Using  the  angular 
representation  [12]  for  r,  [11]  may  easily  be  written  as: 

fc  +-r.f\ 

XT  cos  B  +  YT  sin  B  -  f  - - - 2-  =  -C  [13] 

\  "  J 

By  so  writing  [11],  we  see  that  each  measurement  of  B  and  f  ,  together 
with  the  known  values  of  C  and  Vg,  determines  a  set  of  coefficients 

c  +  f  •  % 

cos  B,  sin  B,  -  - - - C 

a 


for  a  linear  eauation  in  the  three  unknowns  X^,  Y^,  f. 

Thus,  taking  bearing  and  freauency  measurements  simultaneously  at 
N  different  moments  results  in  a  system  of  N  simultaneous  linear 
equations  in  the  three  unknowns. 

Three  such  eauations  would  generally  suffice  to  determine  values 
for  X^,  Y^,  f,  and  the  remaining  N-3  equations  would  either  comDlement 
the  solution  with  redundant  information,  or  else  be  inconsistent.  Note 
that  in  seeking  fixed  solutions  XT»  YT  and  f,  we  implicitly  assume  that 
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during  the  measurement  period,  the  source  changes  neither  velocity  nor 
frequency. 


The  usual  situation,  measurements  with  error  components  (noise), 
will  normally  give  rise  to  an  overdetermined  or  inconsistent  system  of 
equations  of  the  form  C 1 3 3 •  It  is  possible,  in  this  case  to  use  the 
generalized  inverse  of  the  matrix  M,  below,  to  obtain  the  least  squares 
solution  to  the  system: 

•  • 

XT  cos  Bx  +  Yt  sin  B1 

XT  cos  B2  +  Yt  sin  Bg 


•  • 
xt  cos  bn  +  yt  Bin  B: 

or 


A 

fA 

• 

, ?T, 

w 

[lU'l 


where  [ 1 4 ' ]  is  the  matrix  equation  counterpart  of  the  system  [14]. 
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The  theory  of  the  eeneralized  inverse  can  be  found  in  Ref.  1. 
For  example,  in  the  APL  language  used  for  the  programs  listed  in  this 
report,  the  operation  Islh  always  yields  the  generalized  inverse  of  M, 
unless  M  has  fewer  than  3  (in  this  case)  independent  columns.  Thus, 
where  EG  is  used,  the  eeneralized  inverse  always  results;  it  is  identical 
with  the  true  inverse  where  one  exists. 

Solving  the  system  [14]  throueh  the  eeneralized  inverse  as  in  the 
first  part  of  the  APL  function  D0PPVEC2  (Appendix  A)  yields  an 
unconstrained  least  squares  solution  for  source  velocity  and  frequency. 
However,  the  search  problem  is  seldom  unconstrained.  Several  types  of 
constrained  problems  and  their  solutions  are  discussed  in  the  next 
section. 

3.1  Constraints 

The  first  type  of  constraint  provided  by  additional  (outside) 
information  is  easily  disposed  of.  If  frequency  is  known,  velocity 
components  are  found  by  solving  system  [14],  as  we  discussed  above,  and 
this  operation  is  simplified  because  there  are  only  two  unknowns  left, 
Instead  of  three.  The  known  freauency  constraint  is  not  described 
further. 

The  second  type  of  constraint,  that  of  known  direction  of  source 
movement,  without  considering  speed,  is  also  easy  to  handle.  When  the 
direction  is  known,  YT  becomes  a  constant  multiple  of  XT,  and  can  be 
replaced  in  system  [14],  leading  again  to  a  system  to  be  solved  for  two 
unknowns.  If  direction  is  known  to  lie  within  a  certain  angular  sector, 
the  constrained  solution  is  easily  obtained  from  the  unconstrained 
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solution  as  discussed  under  the  next  case. 

The  third  type  of  constraint,  that  of  known  or  bounded  source 
speed,  will  be  carefully  treated  for  its  solution  contains  the  methods 
for  handling  the  other  constraints  mentioned  above. 


3 .2  Constrained Speed 


The  fundamental  linear  system  [14],  which  can  usually  be 
satisfied  only  approximately  in  least  sauares  sense,  can  be  thought  of 
in  the  following  way.  Consider  the  left  side  of  each  eauation  as 
defining  a  coordinate  function,  Z^  ,  of  X^,  Y^,  f,  as  follows: 


+  r  .  VP 


•  • 


D.  (XT,  Yt>  f)  =  ^  cos  Bi  +  Yt  sin  B^fl - j- 


=Z, 


[15] 


If  there  are  N  measurements,  i  varies  from  1  to  N,  and  [15] 
yields  N  coordinate  functions,  which  define  a  mapping  of  the 
three-dimensional  parameter  space  XT»  Y^. ,  f,  into  an  N-dimensional 
space.  The  least  squares  solution  of  system  [14]  can  be  considered  from 
a  geometrical  point  of  view. 

The  least  squares  solution  vector,  <xjj!,  Y®,  f°>  (or  point  in  the 

0  0  0 

three-dimensional  space),  determines  a  vector  <Zy  Z^>  which  is 

as  close  as  possible  to  the  constant  vector  <-C,  -C,  ...,  -C>,  of  N 

0 

coordinates.  If  the  minimum  achievable  value  of  distance  between  <Z ^ , 
...,  Z^>  and  <-C,  ...,  -C>  is  denoted  SQ,  i.e.  the  least  souares 
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0  0 

difference  is  Sn,  then  the  point  <Z  ,  ....  Z  >  lies  on  the  surface  of  an 

IN 

N-airaensional  sphere  with  its  center  at  <-C,  . ..,  -C>. 


0  0 

Point  <Z^ ,  . ..,  Z^>,  at  distance  SQ  from  <-C,  . ..,  -C>  mivht  not 

be  attained  because  of  constraints  on  one  or  more  of  X.^,,  f.  The 

point  at  distance  Sg  is  the  unconstrained  minimum;  if  it  is 
unattainable,  the  algorithm  should  seek  a  solution  in  the  constraint 
region  of  the  three-dimensional  space  of  the  parameters  X^.,  YT,  f  which 
yields  a  value  <Z^,  Z^>  as  close  as  possible  to  <-C,  -C>,  i.e. 

a  constrained  minimum.  Thus,  the  alaorithm  must  seek  a  parameter  DOint 
<xij,,  Y.J,  f%  for  which  the  associated  <Z^,  Z^,  Z*>  lies  on  the 

surface  of  some  other  sphere  about  <-C,  ...,  -C>,  of  radius  S^>Sq.  The 

equation  of  the  sphere  of  radius  is: 


2  ? 

I  {Z.  -  <-C>>  »  S  2 

i  =  1 


UNCLASSIFIED 

16 


Squaring  the  left  side  of  [16]  yields 


A11XT  +  A22YT  +  A33f  +  2A12XTYT  “  2A13XTf 


12  T‘T  13  T 


+  2CA1Xt  -  2A23YTf  +  2CA2Yt  -  2CA3f  +  E  =  0 


where 

N  2  N  2  N  g 

A,  =  E  cos  B.  ,  A  *  l  sin  B.  ,  A„  =  I  E. 

11  i=l  1  i»l  1  i=l 

A12  ®  I  cos  B^  sin  B^ ,  A13  ■  I  ccs  B^  -  Z  sin  Bi 

A^  ~  l  cos  B^,  Ag  =  E  sin  B^  A^  -  I  Pi 

2  2 

E  =  NC  -  S* 

Equation  [17]  represents  a  family  of  ellipsoids  centered  on  the 
•  _  *  . 

point  <X^,Y°,f°>.  The  family  parameter  S^  determines  the  size  of  each 
ellipsoid.  For  S1  =  SQ  there  is  the  point  ellipsoid  at  <X^,Y®,f°>. 
Figure  2  shows  the  geometry  described  above. 


In  the  three-dimensional  space  of  X^,,  Y^,  f,  a  constraint  of 
constant  speed,  V,  is  described  by  the  surface  of  a  cylinder 


'2  ’2  2 

4  +  ■  v 


and  a  constraint  of  bounded  speed  by  the  surface  and  interior  of  this 
cylinder, 


d  +  4  *  v2. 
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For  those  constraints,  the  constrained  minimum  is  obtained  when 
the  solution  ellipsoids  ere  expanded,  by  increasing  S^,  until  tendency 
is  achieved  with  the  constraint  cylinder.  Of  course,  if  the  constraint 
region  consists  of  the  cylinder  and  its  Interior,  and  the  unconstrained 
minimum  X^,  y|,  f®  lies  in  the  cylinder,  then  =■  S^. 

The  mathematical  problem  of  finding  the  point  of  tangency  between 
ellipsoid  and  cylinder  can  be  reduoed  to  two  dimensions,  for  this 
tangency  occurs  exactly  when  tangency  occurs  in  the  X^-Y^  plane  between 
the  circle  X^  ♦  Y2  a  V*  and  the  ellipse  which  is  the  ellipsoid's 
projection  (shadow)  in  that  plane. 

The  projection  of  the  ellipsoid  is  deseribable  as  the  set  of 
points  XT,  Yt  above  which  exactly  one  point  on  the  ellipsoid  can  be 
found.  Equivalently,  if  [173  is  regarded  as  ar  equation  to  be  solved 
for  f,  given  XT,  YT,  the  projection  points  are  those  for  which  the  two 
quadratic  solutions  for  f  coincide.  In  other  words,  they  are  just  the 
points  Xy,  Yj  for  which  tne  discriminant  of  [17],  regarded  as  a 
quadratic  in  f,  vanishes.  Setting  th*t  discriminant  equal  to  zero 
gives: 


(a13xt  -  a23yt  +  GA3)2- 
A33  +  A22YT  +  2Al2*TYT  +  2CA1XT 


2CA2YT 


■f  E)  *  0 


or 


BU*?  *  B22  *T  *  2IWt  +  2BA  *  2B2  *  F  ■  0 


[18] 
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where 

2  2 
B11  *  a13  “  A11  A33  *  B22  ■  A23  ■  A22  A33 

B  ^  *  C  (A13A3  ~  A33A? ) ,  Bg  ■  C  (A23A3  “  A33A2^ 

b12  «  (*13A23  "  A33A12J»  F  C  *3*  A33E 

Equation  [18]  ia  the  ellipse  eouation  to  be  solved  for  tendency 

•  2  *  2  2 

with  the  oirole  equation  XT  +  ^  «  V  .  Finally,  to  ensure  tangenoy, 
another  equation  can  be  added  by  differentiating  the  oirole  equation  and 
[18],  and  then  equating  slopes  at  the  solution  point.  If  this  is  done, 
the  following  system  of  equations  in  three  unknowns  must  be  solved. 
Recall  that  F  is  a  variable  sinoe  it  oontains  the  radius  S1  whioh  must 
be  increased  to  produoe  tangenoy.  S1  appears  nowhere  else. 

+  B22Y^  +  2312XtYt  +  2B1Xt  +  2B2Yt  +  F  -  0 

+  Y2  -  V2  -  0  [19] 

*12*1  -  W?  *  (B22  -  BU>  Vt  *  B2*T  -  BA  '  0 

The  system  [19l  is  solved  by  a  Newton-Raphson  iteration  in  the 

*  1  ’1 

APL  function  D0PPVEC2.  When  the  XT,  Y^. ,  S1  coordinates  of  the  DOint  of 
tangency  are  found,  the  value  of  f  is  recalculated  usina  the  ellipsoid 
[1?].  The  result  of  the  calculation  is  stored  in  an  APL  vector  called 
STV,  whioh  consists  of  f^,  X.J,  . 
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This  completes  the  description  for  a  cylindrical  speed 
constraint . 

3.3  Constrained  Direction 

To  solve  the  constrained  direction  problem  discussed  earlier,  we 
must  only  reconsider  the  geometry  of  the  constrained  minimum  problem,  as 
in  Fig.  3.  In  those  terms,  a  direction  constraint  region  consists  of  a 
(double)  wedge  with  interior,  parallel  to  the  f-axis.  The  boundary 
planes  of  the  constraint  region  are  defined  by  limiting  direction  lines. 

Again,  if  the  unconstrained  minimum  solution  <X^,  Y^,  f°>  lies 
within  the  constraint  wedge,  the  constrained  problem  is  solved.  If  not, 
the  situation  is  identical  to  the  one  discussed  above,  and  S1  must  be 
increased  until  an  ellipsoid  becomes  tangent  to  the  constraint  wedge. 
However,  it  is  not  necessary  to  repeat  the  process  used  for  the  cylinder 
constraint  region.  It  will  be  noted  that  the  desired  point  of  tangency 
is  obtained  when  an  ellipsoid  of  the  family  is  tangent  to  one  of  the 
wedge  boundaries.  Thus  the  constrained  solution  is  one  of  the  two  for 
which 

Yt  *  K^t  or  *T  s  K2Xt  [20] 

It  was  noted  earlier  that  when  a  precise  relation,  like  [20],  is 
obtained  between  Y^,  and  X^,  the  fundamental  system  tl4]  reverts  to  a 
system  of  equations  in  only  two  unknowns,  and  the  least  sauares  solution 
can  be  obtained  using  a  simplified  generalized  inverse. 
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Thus  in  the  oase  of  (wedge)  direction  constraint,  the  fundamental 

*  k 

system  [14]  can  be  solved  twice,  once  with  YT  s  X,p,  and  once  again 
with  ¥t  =  K2xt;  the  best  of  the  two  solutions  (least  squares  sense)  is 
the  solution  of  the  constrained  direction  minimization  problem.  Of 
course,  if  the  direction  constraint  wedge  is  a  plane,  only  one  suoh 
solution  results. 

3.4  Constrained  Direction  and  Velocity 

The  most  realistic  constraint  case  is  perhaps  one  in  which  the 
source  direction  is  known  to  lie  between  two  intersecting  half-lines  and 
the  source  speed  between  two  values  In  this  case,  the 
constraint  region  is  a  cylinder  parallel  to  the  f-axis,  as  in  Figs.  2 

%  k 

and  3,  but  having  for  its  base,  in  the  X^-Y^  plane,  a  figure  oounded  by 
four  ares  which  arc  the  intersections  of  an  annular  ring  with  an  angular 
sector,  measured  from  the  center  of  the  ring. 

Again,  if  the  unconstrained  solution  lies  within  the  constraint 
region,  the  problem  is  solved.  If  it  does  not,  the  best  constrained 
solution  will  then  lie  on  the  boundary  of  the  constraint  region,  at  the 
point  where  the  smallest  ellipsoid  described  by  eauation  [17]  becomes 
tangent  to  the  constraint  region,  by  increasing 

Although  finding  the  point(s)  of  taneeney  in  this  ease  reauires 
more  care  than  it  did  in  the  simple  circular  case  of  speed  constraint, 
it  is  not  difficult  to  develop  an  algorithm  which  will  do  so.  The  fact 
that  the  ellipsoids  of  the  family  described  by  [17]  never  change  their 
center  or  their  principal  axes  as  S1  is  increased  helps. 


at  tfr 
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4.0  ESTIMATION  OF  ftAtGE 

4 . 1  Rough  Ranee  Estimation 

The  first  method  discussed  depends  upon  certain  approximations  to 
the  ranee  vector  in  the  fundamental  eouation  for  apparent  freouencv 

f  (C  +T.7j  =  f  (C  +-J-.-VJ  [21] 

a  i  c> 


The  approximations  will  be  better  for  short  tracking  intervals  and  lone 
ranges.  The  principal  limitation  of  the  rough  ranee  estimate  is  error 
in  frequency  estimation,  f.  With  an  external  input  of  freouencv,  the 
method  could  be  auite  accurate  but,  if  frequency  is  badly  determined, 
the  rough  range  estimate  will  be  useless.  Further,  as  discussed  below, 
the  roueh  range  estimate  is  unreliable  when  both  source  and  searcher 
travel  alone  or  near  the  line  which  joins  them. 


Eouation  [21]  can  be  rewritten  as 


and  f  /f  may,  quite  acceptably,  be  approximated  by  1,  leading  to: 

Cl 


[22] 


1 1 
W 


f 
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To  simplify  the  discussion,  we  will  assume  that  Vg.  is  constant  over  the 
measurement  period  or,  eauivalently ,  that  a  subset  of  the  measurements 
taken  when  V_  is  constant  is  used  for  the  estimate.  A  slishtly  more 
elaborate  procedure  can  be  developed  in  which  Vs  can  vary  and  all 
observations  can  be  used . 

Summing  [22]  over  the  N  usable  observations  yields 

fa  -  f 

N  3i  ^ 

C  l  - 7 -  =  (Vs  •'  VT)  .  S  r.  [23] 

i=1  1  i=1 


This  summation  reduces  the  importance  of  independent  random  errors  in 

the  measurements  of  the  f  ,  but  an  examination  of  the  effect  of  the 

ai 

initial  range,  R^,  on  the  rieht  side  of  [23]  reveals  the  sensitivity  of 
this  equation  to  determination  of  the  freouency  f. 

If  R i  is  larse  and  the  time  interval  between  measurements  is 
relatively  small,  there  will  be  little  chanve  in  the  unit  direction 
vectors  r^  ,  and  we  can  write 

N  __ 

£  r  =  N  r  [24] 

i=  1 

where  r1  is  the  initial  unit  direction  vector.  On  the  other  hand,  if  R1 
is  close  to  0,  the  sum  of  the  unit  vectors  is 

*v’  . 

l  r±  =  v  +  (H  -  1)  _ T _ s_ 

1=1  I  I  \  %  I  I 


[25] 
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for  the  first  measurement  is  taken  in  the  direction  r1  of  the  source, 
but  all  subseauent  ones  are  in  the  direction  of  the  relative  velocity 
V^-Vg.  (More  complex  expressions  are  clearly  needed  here,  and  below,  if 
Vg  is  allowed  to  vary).  The  expressions  [24]  and  [25]  can  be  used  to 
bound  values  on  the  rieht  side  of  [23],  as  follows: 


"r’-j  •  (Vs  -  VT)  -  (N  -  1)  |  | V  -  VT|  |  S  (Vg  -  VT).  £  PJ_ 


N 

£ 

1  =  1 


5  "  '•r  <’s  -  V 


[26] 


where  ||V  “  VT | (  is  the  magnitude,  or  IchikLIi  of  the  vector  V<,  -  V^,  and 

is  positive.  It  is  clear  that  the  lower  bound  in  [26]  is  always 
inferior  to  the  upper  one,  although  both  may  be  negative .  The  upper 
bound  is  negative  onlv  when  the  (signed)  projection  of  the  vector  V^-V^ 
onto  the  direction  r^  is  negative,  or  when  the  (siened)  projection  of 
VT  -  Vg  onto  r^  is  positive,  or  when  the  source  is  moving  away  from  the 
searcher . 


Since  [23]  can  be  used  to  replace  the  middle  expression  in  the 
inequality  [26],  the  following  inequality  results: 


(N  -  1) 


11%  -%l  I  *  C 


N 

£ 

i  =  1 


5  N  r 


1 


[27] 


It  is  in  the  form  [27]  that  the  sensitivity  of  the  roueh  ranve 
estimate  (still  to  be  described)  is  most  evident.  If,  for  example  V 

^  i  O 

lies  close  to  the  line  of  action  of  r^,  then  r^  .  (Vg  -  V^.)  is  neeative, 
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and  close  in  value  to  | |V_  -  V_ 1 1  .  In  this  case  the  two  bounds  in  [27] 

o  T 

differ  by  little.  Then,  if  the  estimate  of  VT  is  trustworthy,  [27] 
might  be  used  to  provide  a  worthwhile  estimate  of  f;  however,  unless  f 
is  well  known,  [27]  will  not  be  useful  in  finding  the  value  of  R^. 

If  the  bounds  in  [27]  are  spaced  far  enough  apart,  and  if  a 
reliable  estimate  of  f  is  available,  a  method  for  routrh  ranee 
determination  results  from  closer  examination  of  the  Quantity 

—  N 

(vs  "  V  •  T-  =  Q(Ri}  C28:i 


The  inequality  [?6]  shows  the  maximum  possible  variation  of  0(R.)  as  R^ 
varies  from  0  to  «  and,  if  the  measured  values  f  satisfy  the 
ineauality  [27],  a  value  of  R,  can  be  found  to  make  true .  Such  a 

value  of  R.j  would  be  the  rough  range  estimate. 


To  examine  the  behavior  of  0(R.j)  as  r1  is  varied,  two  other 
approximations  are  made  to  the  unit  vectors "r 1 .  If  TT^  represents  the 
time  at  which  measurement  i  ii  made  after  the  first  measurement,  i.e. 
TT1  =0,  then 


\  +  TT±  (Vt  -  V 
\\\  +  TTi  (VT  -  %)\\ 


[29] 


It  should  be  noted  here  that  the  time  TT^  is  not  modified  to  obtain  the 
position  of  the  source  at  the  moment  of  the  emission  of  pulse  i  rather 
than  at  the  time  of  its  reception.  Such  corrections  are  made  later  in 
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this  section  as  well  as  in  computer  programs  to  allow  for  the  source 
movement  during  transit  time  of  the  pulse.  Here,  however,  all  times  are 
measured  from  the  time  of  reception  of  the  first  pulse  so  that  onlv  a 
minimal  correction  would  be  needed  for  the  small  change  of  range  between 
the  first  measurement  and  measurement  i . 


Further,  approximating  the  denominator  in  [29]  by 
♦  OT1‘?1  •  <VT  -%) 


which  represents  the  value  of  corrected  by  the  change  in  range 
produced  by  the  projection  of  relative  velocity  in  the  r1  direction, 
results  in  an  approximation  to  [29]: 


+  TTi  (VT  -  Vs) 
+  TT.^. fVT  -  Vg) 


[30] 


Now,  if  [30]  is  substituted  into  [28]  and  if  the  result  is 
differentiated  with  respect  to  R^»  the  following  derivative  is  obtained: 


N 


(?1.(vs-Tt))2)  E 


TTj 


i-1  (ai+TT1r1.(VT-Vs))‘ 


[31] 


and  the  auantity  is  always  Dositive. 

Finally,  when  the  inequality  [273  is  satisfied  by  the  estimated  f 

and  the  measured  f  ,  the  Quantity  0(R  )  rises  steadily  with  increasing 
ai  1 

R^,  one  unambiguous  value  of  R^  can  be  found  for  which  [23]  is 
satisfied.  This  value  is  the  rough  range  estimate.  The  easiest  wav  to 
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find  it  is  to  apply  the  method  of  Newton-Raphson  to  [30.  The  estimates 
are  found  with  the  APL  function  RRNG  of  Appendix  A. 

One  further  look  at  the  form  of  the  derivative,  in  [31],  reveals 
again  the  sensitivity  of  this  estimate  to  the  direction  of  r^  and  the 
relative  velocity  VT  -  Vg.  Since  the  derivative  contains  the  factor 

(11%  -  vT||2  -  %  .  (Vs  -  v“T))2) 

it  may  be  very  small  when  VT  -  Vg  lies  in  the  line  of  action  of"^  •  The 
range  estimate  then  becomes  auite  unreliable,  and  the  Newton-Raphsou 
root- finding  techniaue  may  have  to  be  abandoned  for  another  teehniaue 
which  is  described  below. 


4.2  Best  Linear -Range  Hatsii 

The  second  method  for  possible  ranne  determination  also  uses  the 
least  squares  solution 


STV  = 


-  <f.  v  v 


[32] 


of  the  linear  system: 


c  +  ri  •  V 


XT  cos  B1  +  Yt  sin  B1  -  f  - — 


=  -C 


C  +  rN  *  V 


XT  008  bn  +  yt  Bin  Bn  -  f  — F 


=  -c 


[33] 
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In  finding  the  least  souares  solution  STV,  the  actual  measured 
values  of  bearing  B,,  r  were  used  to  evaluate  the  coefficients  of 

system  [33].  Once  values  of  f,  X^,,  'ij  are  available,  one  can 

a)  assume  a  value  of  range  at  initial  measurement,  R-,; 

b)  use  the  measured  initial  bearing,  the  estimated  velocity  XT, 

together  with  the  assumed  value  of  to  venerate  a  set 

of  bearings  B_,  B*,..,  B  at  the  later  measurement  times; 
c.  5  N 

c)  calculate  the  left  sides  of  the  eauations  in  [33]  and  compute 

the  summed,  squared,  differences  between  left  and  right  sides 

in  [33],  using  the  results  of  b)  and  the  other  known  and 

measured  quantities  C,  V„  ,  f  . 

si  ai 

The  result  obtained  in  c),  called  DELTASO  in  the  APL  functions  of 
Appendix  A,  can  be  considered  a  measure  of  merit  of  the  assumed  ranire 
R.j ,  subject  to  the  accuracy  of  the  estimate  STV.  If  there  is  a  minimum 
value  of  DELTASO  over  the  interval  of  possible  assumed  ranges,  R^  the 
value  of  R^  correspond  ins:  to  the  minimum  could  be  chosen  as  the  best 
linear  range  match. 

The  above  computational  procedure  is  easily  implemented.  The 
calculation  is  performed  with  the  APL  function  RNSEL  of  Appendix  A. 

4.3  usaaL  amaiiaa ItscatjLan. 

The  third  method  for  range  estimation  is  the  classical  least 
squares  iteration  that  permits  to  find  the  point  in  parameter  sDace, 
<f,  XT,  R^,  which  minimizes  a  sauared-summed  error  function. 
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Like  the  best  linear  range  method  or  second  method  previously 
explained,  the  third  method  uses  values  of  the  parameters  to  generate  a 
set  of  bearings,  so  that  its  comparison  with  measured  data  involves  only 
frequency  measurements. 

In  the  second  method,  only  the  ranae  parameter  can  be  varied.  In 
the  third  method,  the  minimum  is  sought  by  evaluating  the  zero  points  of 
partial  derivatives  (of  the  sauared-summed  error  function).  It  is  thus 
possible  to  vary  any  of  the  parameters. 
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Figure  L  -  Example  of  range  estimation  for  two  searchers 
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Figure  5  -  Example  of  range  estimation  for  one  searcher 

Unfortunately,  in  the  presence  of  measurement  noise,  the 

squared-summed  error  surface  for  freauencies  seems  to  possess  many 

distinct  relative  minima.  Computation  reveals  that  these  are  not  always 

located  near  the  true  values  of  the  parameters  but,  in  the  search  for 

the  zeros  of  partial  derivatives,  incorrect  minima  often  arise.  Fieure 

4  shows  the  result  of  a  typical  calculation  in  which  a  search  was 

performed  on  a  range  arid,  and  the  algorithm  was  allowed  to  seek  the 

*  * 

best  values  of  f,  X^,,  1^,  to  minimize  SIGMASO  in  the  APL  function 

KITERATE.  In  this  case,  two  searchers  are  used  and  the  ranee  is  not 
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exoessive;  so  it  is  a  favorable  case.  The  figure  shows  that,  without  a 
good  initial  range  estimation,  an  algorithm  searching  in  range  as  well 
as  in  the  other  parameters  might  yield  several  different  values.  The 
behavior  pictured  is  not  exceptional. 

The  case  presented  in  Fig.  5  was  calculated  for  one  searoher. 
Although  the  single  searoher  was  allowed  to  change  course,  and  the  same 
total  number  of  observations  was  used  as  in  the  two-searcher  case  of 
Fig.  4,  it  is  apparent  that  many  more  local  minima  exist.  Countine  end 
point  minima,  there  are  at  least  five  in  Fig.  5,  and  an  iteration  might 
converge  to  any  one. 

Despite  the  complexity  of  the  function  to  be  minimized,  the  least 
squares  algorithm  is  useful  in  many  cases  and  it  often  yields  a  deepest 
minimum  not  far  from  the  true  vslue  of  the  ranee  parameter.  The  basic 
equations  for  the  method  are  presented  next. 

Using  [21],  solved  for  the  apparent  freaueney,  f  ,  the 
squared-summed  differences  for  the  frequency  observations  can  be 
expressed  as 


N  C  +"rt 

SIGMASQ  *  E  (f  - - 

i=l  C  +  r* 


) 


2 


[34] 
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SIQMASQ  is  to  b«  minimised  by  fixing  a  range,  then  seeking  points  where 
the  partial  derivatives  below  vanish. 


.  »  C  r,  .  Ve  C  ♦  r,  .  V- 

-  2  e  (f  — -r  ■  -  f  )  — —■  ■  J 

3f  i-1  C  +  p1  .  VT  ai  C  +  r,  .  Vm 


i  C  +  T±  .  VT 


K 

ir-  •  2  Z  ( 
i-1 


3^ 


[35] 


N 

2  Z  ( 


i-1 


where  the  fa,  on  the  right,  is  the  function  of  XT  and  YT  obtained  by 
solving  [21]  for  f#.  However,  these  formulas  must  be  interpreted 
carefully.  For  each  i,  1<  1  <  N,  fa(XT»  is  a  different  function 
and,  hence,  each  equation  of  [35]  is  a  scaler  product  between  data 
vectors,  of  length  N.  It  is  useful  to  recognize  this  fact  in  notation, 
by  using  a  scaler  product  notation  between  data  vectors  of  length  N, 
and  by  emphasizing  such  vectors  by  the  notation  — Then,  if  the  error 
vector  whose  components  appear  in  the  parenthesis  on  the  right  sides  in 
[35]  is  denoted  EV,  [35]  can  be  rewritten  in  compact  form  as  follows: 
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3f 

Then,  —A  is  an  N-veotor  whose  ith  component  is: 
3XT 


(C+ri-VT)‘ 


5«(N?i-VT)(V8  --rtMc^-V  )(QM  [36] 

i  OJVp  1  OArp 


3^ 
9f 


and  there  is  a  similar  expression  for  the  ith  component  of  —Ji.  . 

3YT 

Thus,  the  evaluation  of  the  needed  partial  derivatives  oannot  be 
completed  without  evaluating  the  expressions 


d  r. 


3xt 


and 


Sr. 


3Y„ 


It  must  now  be  recalled  that  r^  is  the  unit  vector  pointing  from  the 
position  of  the  searcher  at  time  i,  TT^  *  toward  the  position  of  the 
source  at  a  time  TTR^  at  which  the  ith  pulse  was  emitted.  In  other 
words,  the  source  position  times  must  be  corrected  for  transmission 
time.  This  correction  is  performed  with  the  APL  function  NUVECT,  for 
use  in  the  iteration  function  KITERATE. 


In  terms  of  TT^  and  TTR^,  r^  can  be  written  as: 


<R1 

cos 

B1 

+  XT 

TTRi 

R, 

sin 

B, 

+  Y_ 

TTR, 

1 

1 

T 

i 

Xs  TTi, 
* 

Y_  TT.> 

b  1 


[37] 


where  R^  is  the  magnitude  of  the  vector  in  corners  <>. 
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Next,  [37]  is  differentiated  with  respeot  to  XT  and  YT,  to  yield 

the  two  needed  derivatives.  Note  that  no  approximation  has  been  made  in 

obtaining  these  derivatives: 

3r.  TTR.  _  2 

-  - *•  <RRN*  -  RRX  ,  -  RRX  x  RRY> 

ax,,  rrn-5 

[38] 

3r.  TTR^  2  2 

——  »  - i  <-RRX  x  RRY ,  RRN  -  RRY  > 

3Yt  RRNJ 

where 


a  • 

RRX  ■  R.^  cos  B1  +  TTRj^  x  X^  -  x  Xg 

•  e 

RRY  -  R^^  sin  +  TTR^  x  Yj  -  TTj  x  Yg  L  ->9  ] 

RRN  «  (RRX2  +  RRY2)^ 


and  these  Quantities  represent  the  X  and  Y  coordinates  of  the  source  as 
well  as  distance  in  a  corrected  searcher  relative  space  in  which  the 
source  position  at  the  time  of  emission  of  pulse  i  is  related  to  the 
searcher  position  at  the  time  of  reception  of  pulse  i. 


The  expressions  [38]  substituted  into  [35],  substituted  into 
[35']  yield  formulas  for  the  needed  partial  derivatives.  It  is 
pointless  to  write  out  these  complete  formulas  here.  The  reauired 
quantities  are  computed  in  the  APL  function  FCN. 


Finally,  the  search  for  the  zero  DOints  of  the  derivatives  of 
[35']  is  carried  out  by  a  Newton-Raphson  procedure  for  which  the 
required  derivative  is  calculated  numerically  with  the  APL  function 
JACOB,  and  the  actual  iteration  to  the  convergence  criterion  is 
performed  with  KITERATE.  The  search  is  managed  by  a  control  function 
named  SSSTM. 
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5.0  USE  OF  ALGORITHMS 

Computer  programs  are  available  for  only  two  of  the  algorithms 

•  * 

analyzed  above.  DOPPVEC2  estimates  f,  XT  and  with  speed  oonstraint. 
SSSTM  estimates  target  range  by  least  souares  iteration  from  the 
D0PPVEC2  results  and  then  makes  a  new  estimation  of  f,  XT  and  Y^.  These 
computer  programs  are  listed  in  Apprendix  A. 

If  the  reader  uses  these  programs,  he  may  notice  a  few 
peculiarities  of  behavior  which  might  usefully  be  pointed  out  here: 

1)  If  bearing  and  freouenoy  measurements  are  totally  free  of 
error,  DOPPVEC2  will  yield  an  exact  estimation  of  f,  X^,  Y^..  SSSTM, 
however,  will  yield  estimates  with  small  errors  because  the  search  in 
range  is  carried  out  at  discrete  *alues  which  may  not  correspond  to  true 
range. 


2)  The  variable  SIGMASQ  does  not  always  pass  through  minima;  it 
may  decrease  indefinitely  as  range  increases  in  the  interval  searched  by 
SSSTM.  In  such  cases,  the  estimated  value  of  range  is  always  near  the 
high  end  of  the  interval  and  mav  bear  no  resemblance  at  all  to  the  true 
value . 


3)  In  principle,  SSSTM  improves  the  estimate  of  f,  XT»  and  Y^ 
made  by  D0PPVEC2.  In  reality,  the  new  estimate  may  be  worse  than  the 
original . 

4)  Both  algorithms  favor  speed  estimates  close  to  the  speed 
constraint  rather  than  to  the  real  speed.  In  most  cases,  the 
constrained  solution  is  a  bounding  solution  in  the  constraint  region. 


The  operational  effectiveness  of  the  algorithms  discussed  above 
is  unknown.  It  can  be  estimated  by  simulation  studies  in  which  searcher 
tactics  are  taken  into  account.  Such  studies  have  now  been  undertaken 
for  the  two  algorithms  already  programmed.  The  other  algorithms 
analyzed  above  should  be  easy  to  program  whenever  needed. 


UNCLASSIFIED 

38 


Ben-Israel 

Applications 


7.0  REFERENCES 

and  Greville,  "Generalized  Inverses,  Theory  and 
,  New  York,  Wiley-Interscience,  1971*. 


UNCLASSIFIED 

39 

APPENDIX  A 

THE  APL  FUNCTIONS 

The  following  APL  functions  are  presented,  and  briefly  described. 

1)  MEASGEN:  Accepts  the  basic  descriptions  of  source  and 

searcher(s),  and  generates  bearing  and  freauency  inputs,  with 
random  errors,  for  use  by  the  processing  algorithms. 

2)  D0PPVEC2:  Performs  the  least-squares  solution  of  the 

fundamental  system  of  linear  equations,  subject  to  a  speed 
constraint . 

3)  NUVECT:  Uses  the  measurement  times  and  the  first  measured 

bearing  to  convert  an  estimate,  STV,  of  f,  XT,  YT,  and  an 
estimate  of  range,  RN,  into  predicted  source  coordinates  and 
frequencies,  for  comparison  with  the  measurements  generated 
by  MEASGEN. 

M)  RMNL:  Chooses  a  vector  of  N  random  numbers  from  a  Gaussian 

distribution  of  mean  P[1]  and  standard  deviation  P[2].  The 
choices  are  independent. 

5)  SIGM:  Computes  the  summed-squared  deviation  for  freauency 

data,  SIGMASQ. 

6)  FCN:  Computes  the  values  of  the  partial  derivatives  of 
SIGMASQ,  with  respect  to  the  parameters  f,  X^.,  Y^,. 
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7)  JACOB:  Performs  a  numerical  calculation  of  the  matrix  of 

second  partial  derivatives  of  SIGMASQ,  for  use  in  the 
iteration  KITERATE. 

8)  KITERATE:  Does  the  Newton-Raphson  iteration  to  find  a 

zero-point  for  the  partial  derivatives  of  SIGMASQ. 

9)  RRNG:  Computes  the  rough  ranee  estimate. 

10)  RNSEL:  Chooses  the  best  linear  ranee  estimate  for  the 
solution  of  the  fundamental  linear  system. 

11)  SSSTM:  Performs  a  grid  search,  using  KITERATE,  for  a  ranee 

at  which  the  smallest  value  of  SIGMASQ  can  be  obtained,  by 
simultaneously  varying  the  other  parameters. 

12)  DAPPROX,  DEXACT,  SHIFT  and  SHIFTY:  Solve  the  exact  and 

approximate  formulas  for  the  Doppler  shift,  as  presented  in 
Section  2.0. 

Further  comments  accompany  certain  of  the  individual  APL  functions. 

Function  1 :  MEASGEN 

INPUTS:  TAR  is  a  1  x  8  vector  which  consists  of  the  following: 

TAR[ 1 ]  Source  initial  X 

TAR[2]  Source  initial  Y 

TAR[3]  Source  X  velocity 

TAR [ 4 ]  Source  Y  velocity 
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TAR [5  3  Maximum  possible  source  speed 

TA R [ 6 ]  Default  estimate  of  source  range  (to  be  used  if  all 
efforts  to  find  range  fail) 

TA R C 7  3  Source  emission  frequency 
TAR[8]  Speed  of  sound  in  the  medium 

Since  only  one  speed  vector  is  given  for  it,  the  source  moves  on 
a  straight  line,  at  least  throughout  the  measurement  interval. 


SRCH1  is  a 
SRCH 1 [ 1  ] 
SRCH 1 [2] 
SRCH 1 [3] 
SRCH 1 [4] 
SRCH 1 [5] 
SRCH 1 [6] 
SRCH 1 [7] 


1  x  7  vector  which  consists  of  : 
Searcher  initial  X 
Searcher  initial  Y 
Searcher  initial  X  velocity 
Searcher  initial  Y  velocity 
Searcher  turn  time 
Searcher  final  X  velocity 
Searcher  final  Y  velocity 


SRCH2  is  a  1  x  7  vector  with  the  same  components  as  SRCH1, 
describing  a  second  searcher.  Although  only  two  searchers  are  provided 


for  and 

only 

a  simple  dogleg  path  is 

given 

for  each 

one 

there  are 

no 

limits 

in 

the  algorithms  which 

would 

prevent 

reprogramming 

for 

additional  searchers,  and  more  complicated 

search  paths. 

Searchers 

can 

be  prevented  from  maneuvering  by  setting  SRCH1  [5]  and  SFCH2  [5]  longer 
than  the  last  measurement  time  in  TIMEM,  and  the  present  programs  can  be 
used  for  either  1  or  2  searchers,  as  controlled  by  the  input  TIMEM, 
detailed  below. 
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MEASM  is  a  2  x  5  matrix  of  measurement  parameters,  as  follows: 


MEASM[1 ; 1  ] 

MEASM [2;  1] 

MEASM [1 ;2] 

MEASM [1 ; 3] 

MEASM [1 ;4] 

MEASM [1 ;5] 
MEASM[2;2]-[2;5] 


1 

2 

The  mean  of  the  bearing  error  distribution 
The  standard  deviation  of  bearing  error 
The  mean  of  the  freauency  error  distribution 
The  standard  deviation  of  frequency  error 
Same  as  MEASM  [1;2]  -  [ 1 ; 5 ]  for  the  second 
searcher;  these  parameters  can,  of  course, 
be  given  different  values  for  each  searcher. 


TIMEM  is  an  N  x  2  matrix  which  specifies  the  times  of 
measurements.  The  first  column  specifies  the  number  of  the  searcher 
that  makes  the  measurement,  and  the  second  one  gives  the  times.  While 
the  present  algorithms  do  not  consider  loss  of  contact,  this  factor 
could  be  introduced  by  producing  a  random  generator  for  contact  times 
and  searcher  numbers,  or  by  accepting  TIMEM  as  an  output  from  a  search 
algorithm . 

To  use  MEASGEN  for  a  single  searcher,  the  first  column  can  be  set 
entirely  to  1  (or  to  2),  or  else  vectors  SFCH1  and  SFtCH2  can  be  made 
identical . 


FUNCTION:  Lines  [1]  through  [11]  set  the  measurement  times, 
together  with  the  searcher  coordinates  and  velocities  at  those  times. 
True  measurement  times  are  used.  Lines  [12]— [153  find  source  relative 
coordinates  in  searcher  relative  space(s)  at  measurement  times.  Lines 
[16]— [25 ]  perform  corrections  to  obtain  the  source  positions  at  the  time 
of  pulse  emissions,  relative  to  searcher  positions  at  times  of 
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reception.  The  remainina  lines  calculate  true  bearinas  and  freouenoies 
and  lines  [35],  [36]  apply  normally  distributed  independent  errors. 


OUTPUT:  The  output  of  MEASGEN  is  an  N  x  12  matrix  called  DATAM. 
Each  line  stores  the  information  aporooriate  to  a  measurement,  with 
components  as  follows: 


DATAM [I 
DATAM [I 
DATAM [I 
DATAM [I 
DATAM [I 
DATAM [I 
DATAM [I 
DATAM[I 
DATAM [I 
DATAM[I 
DATAM [I 
DATAM [I 


1] 

21 

3] 

4] 

5] 

6] 

7] 

8] 

9] 

10] 
11] 
12] 


Line  number  I 

Searcher  number  for  line  I 
Measurement  time 
Bearine 
Freaueney 

Searcher  X  coordinate  at  measurement  time 

Searcher  Y  coordinate  at  measurement  time 

Searcher  X  velocity,  this  time 

Searcher  Y  velocity,  this  time 

Ranee  from  searcher  1  to  source,  this  time 

Ranee  from  searcher  2  to  source,  this  time 

Turn  time  of  the  searcher  on  this  line 


Function  2:  D0PPVEC2 

INPUT:  DATAM,  from  MEASGEN 

FUNCTION:  Lines  [1]  -  [10]  comDute  the  matrix,  M,  of  the 

fundamental  linear  system  and  obtain  the  least  souares  solution,  usin^ 
the  Moore-Penrose  inverse  ffl  nroerammed  into  APL.  Line  [16]  and  [I'7] 
check  that  the  speed  constraint  specified  in  TAR  [8]  is  met.  If  not, 
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lines  [18]  through  [32]  set  the  geometrical  parameters  of  the  ellipses 


described  in 

the 

text, 

and 

the 

remainder  of 

the 

function 

uses 

Newton-Raphson 

to 

solve 

for 

the 

smallest  value 

of 

DELTASO 

which 

satisfies  the  speed  constraint. 

OUTPUT:  STV  is  the  least  squares  solution  of  the  linear  system, 

subject  to  the  speed  constraint. 

Function  3:  NUVECT 

INPUTS:  DATAM,  from  MEASGEN 

STV,  from  D0PPVEC2,  or  from  KITERATE  when  used 
for  iteration 

RN,  from  SSTM,  whicli  controls  the  grid  search 
for  range. 

FUNCTION:  Lines  [1]  through  [12]  use  the  searcher  information 
Irom  each  line  in  DATAM  to  generate  relative  source  coordinates  and 
distance,  in  searcher  relative  space,  at  the  time  of  reception  of  a 
pulse.  Lines  [133  and  [ 1 4 ]  correct  the  source  position  to  the  emission 
time  of  the  pulse.  After  position  correction,  the  remaining  lines 
compute  the  angular  and  frequency  data  to  be  expected .  The  entire 
operation  is  performed  simultaneously  for  the  whole  of  the  matrix  DATAM. 

OUTPUT:  The  relative  position  and  distance  vectors  (RRX,  RRY, 
RRN),  the  bearing  data  (COSV,  SINV),  and  the  components  of  the  Doppler 
formula  3,  CRVT)  are  all  used  in  other  functions.  Vectors  NUAV,  the 
predicts  observed  frequencies,  and  TTVR,  the  corrected  emission  times, 
are  perhaps  the  principal  outputs  of  NUVECT. 
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Function  6:  FCN 

INPUTS:  DA TAM  plus  the  position  and  freauency  predictions 
of  NUVECT . 

FUNCTION:  Lines  C  3  3 »  [7],  [8]  compute  the  three  partial 

derivatives  of  SIGMASQ,  with  respect  to  f,  XT»  YT,  respectively. 

OUTPUT:  FV  is  a  3  x  1  matrix  of  the  values  of  the  partial 
derivatives. 

Function  7s  JACOB 

INPUTS:  FV,  or  FCN  STV,  from  FCN 

FUNCTION:  Computes  numerical  approximations  for  the  nine  second 

partial  derivatives  of  the  three  first  partials  in  FV,  with  respect  to 

%  , 

each  of  the  three  parameters  f,  XT»  Y^.  Instead  of  a  true  derivative,  a 

-4 

difference  quotient  is  computed  for  an  increment  in  f  of  10  and 
increments  in  X^,,  YT  of  10 

OUTPUT:  The  3  x  3  matrix  JA,  or  JACOB  STV. 

Function  8:  KITERATE 

INPUTS:  DATAM,  TAR,  FCN  STV,  JACOB  STV,  REPIT, 


outputs  of  NUVECT. 
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FUNCTION:  Lines  [4]  and  [5]  compute  lower  and  upper  freouenoy 
bounds  (LNUB,  UNUB)  for  the  frequency  parameter,  based  upon  the  observed 
frequencies  and  the  speed  constraint. 

KITERATE  attempts  to  move  from  point  STV,  in  a  direction 
indicated  by  JACOB  STV,  to  find  another  parameter  vector,  STVA,  which 
satisfies  the  speed  and  frequency  constraints,  and  results  in  a 
decreased  value  of  SIGM  STV.  The  base  values  are  referred  to  as  SIGMB, 
FCNB,  STVB.  If  the  new  value  STVA  of  line  [11]  does  not  satisfy  the 
restrictions  in  [13],  other  values  of  STVA  are  tried  by  progressive 
halving  of  the  parameter  vector  displacement,  [13],  until  one  is  found 
that  satisfies  the  restrictions;  if  it  is  impossible,  th-'n  we  revert  to 
the  previous  estimate,  line  [1?]. 

If  improved  values  of  STV  are  found,  this  process  continues  until 

the  difference  of  summed  squared  deviations,  SIGM  STV,  stabilizes  to 

-5 

within  ±io  until  a  predetermined  number  of  repetitions,  REPIT,  is 
reached . 

OUTPUTS:  STV,  STVA,  both  names  for  the  same  parameter  vector. 

Function  11:  SSSTM. 

INPUTS:  STV,  from  D0PPVEC2. 

FUNCTION:  Lines  [1]  through  [7]  calculate  KITERATE  STV  for  the 

initial  input  STV,  oith  values  of  RN  from  10  to  120.  The  results  are 
stored  in  a  matrix  SSTM,  as  indicated  in  line  [5].  Lines  [8]  through 
[15]  construct  a  vector,  AV,  of  ones  and  zeroes  with  a  one  in  esch 
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position  where  the  corresponding  row  of  SSTM  exhibits  a  relative  minimum 
/alue  of  SIGMA. 

If  there  are  no  relative  minima  in  column  5  of  SSTM,  line  [17] 
prints  a  row  of  zeroes,  as  row  [13]  of  SSTM,  and  terminates. 

Line  [193  selects  the  relative  minimum  rows  of  SSTM  and  places 
them  in  a  matrix  called  RNS.  Lines  [20]  through  [32]  expand  the  search 
about  the  relative  minima  of  SSTM  in  increments  of  -5,  0,  +5  about  each 
relative  minimum. 

These  calculations  use  the  values  of  STV  found  by  KITERATE  at  the 
relative  minimum  ranges,  not  the  original  STV  from  DOPPVEC2.  The 
results  of  the  expanded  calculations  are  recorded  as  new  rows  in  the 
matrix  SSTM. 

Line  [33]  selects  the  absolute  minimum  in  column  5  of  the 
expanded  matrix  SSTM.  Finally,  lines  [34]  through  [45]  further  refine 
the  search  in  increments  of  -2.5,  u,  +2.5  about  the  range  value  found  by 
C 33 3 ,  and  records  the  output  MIN. 

OUTPUT:  MIN. 

Function  12:  DAPPROX,  DEXACT,  SHIFT,  SHIFTY. 

DAPPROX  computes  the  DOPPLER  approximation  of  formula  [7], 
Section  2.0,  as  DAA,  line  [2].  Line  [33  computes  formula  [9'],  as  DAA2. 
DEXACT  solves  the  exact  relation  [5]  for  A  ,  by  a  Newton-Raohson 

Cl 

iteration,  to  accuracy  10”  A  .  The  functions  SHIFT  and  SHIFTY  simDlv 

a 
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permit  the  values  of  exaot  and  approximate  solutions  to  be  computed  and 
displayed  for  a  number  of  different  frequencies,  and  for  a  number  of 
different  velocities  of  source  and  searcher. 

The  vector  input  IV  consists  of: 

Source  X,  Source  Y,  Searcher  X,  Searcher  Y 
N  is  the  desired  number  of  repetitions 
D  is  the  interpulse  interval,  D  =  1/f. 

SHIFT  makes  N  random  choices  of  the  relative  coordinates  X  and  Y, 

SHIFTY  demands  initial  X  and  Y,  as  IV[5],  IV[6],  and  moves  the 
along  a  straight  line  until  N  positions  have  been  achieved. 


while 

source 
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7 NEASG BNlQ 37 

IMEASGENi  TVtLiSlXViSlYVi  S2XV\ S2YVX  RBLXURELY1\RBLX2  iRE 
LI2iBBAR2lRVURV2iFREQ2iNiS 
Cl]  TV^TIMSM  [»2  3 

[2]  TPX+TARlll+TARl3]*TV 

C  3  ]  rpy^JMP^+JviflUSxrv 

[43  SlXP+((SRCRlLl'i+SRCmiy}*TV)*(TV4SRCmlS']))+<.  (SRCHlZll 
+  SRC  H 1  [  3  3  X  SRC  R 1 C  5  3 )  ♦  ( T  V  - SR C  R 1 C  5  ] )  *  SR C  H 1  [  6  3 )  *  ( 2*  V  >  SRC  H 1  [ 
S3) 

[5]  SlYP+((SRCmL2]+SRCHlL*']*TV)x(TVSSRCHll$]) )+( (Soffit 2 3 
+  SR  C  H 1 C  4  3  *  SRC  R 1 C  5  3  )  +  ( T  V  - SRC  H 1 C  5  3 )  *  SR C  H 1  [  7  ] )  x  ( T  V  >  SRC  H 1 C 
S3) 

[6]  S2XP+((SRCH2ll]+SRCH2l3]*TV)*(TV*SRCH2lS]))+( (SRCH2ZH 
+  SRCH2ZSl*SRCH2ZS'})  +  (TV-SRCH2ZS'])*SRCH2ZG))*(.TV>SRCH2Z 
53) 

[7]  S2y?«-(  (SRCH21 2  ]+ SRC  H  21*1*  TV)*  (TV&SRCH21  5 ] ) )+ ( (SflCP2[ 2] 
+  SRC  R  2  C  4  3  x  SRC  R  2  [  5  ]  )  ♦  ( T  V  - SR  C  R  2  C  5  ] )  *  SRC  H  2  [  7  ] )  *  ( T  V  >  SR  C  H  2  C 
53) 

[8]  SlXV+tSRCHlZGlxiTVsSRCinZS])  )+SPCff l[ 6 ]x [ 5 ] ) 

[9]  SlYV‘-(SRCHlZM*(TVSSRCHlZS'}))  +  SRCHlZ7'l*U'V>SRCHlZS )) 

[10  3  S2XV+(SRCH2Z31*<.TVZSRCH2ZS'j))  +  SRCH2ZG]*{TV>SRCH2ZSD 
[11]  S2YV+(SRCH2Z'*]*(TVZSRCH2ZS])  )  +  SRCH2Z7l*[TV>SRCH2ZS']) 

[12  3  RBLX1--TPX-S1XP 

[13]  RBLY1+TPY-S1YP 

[14]  RBLX2+TPX-S2XP 

[15]  RBLY2-TPY-S2YP 

[16]  RV1+((RELXI*2)+RBLY1*2)*.S 

[17]  RV2+-URBLX2*2)  +  RBLY2*2)*.S 

[18]  TPXI+TPX-TARZ31*RV1*TARZB1 
Ci9]  rpy2-*-rpy-r>i/?[3]xpK2>r/iP[8] 

[20]  TPYl+TPY-TARZm*RVHTARZ9 ] 

[21]  rPy2^TPy-Ti4P[4]xpy2*2,i4P[8] 

[22]  PELyi^rpyi-siyp 

[23]  RELY1+TPY1-S1YP 

[24]  RELX2+TPX2-S2XP 

[25]  RELY2+TPY2-S2YP 

[26]  BEARl  +  (  30RELYHRELX1  )*  (RELX\>0  )*  (RELYING  )  )  +  (  (  SORE  l 
Yl*RELXl)+01)x(RELXl <0 ) )+( ( ( " ZORELYl *RELX 1 )+02 ) * {RECX 1 

>o)*(i?PLyi<o)) 

[27]  BEARl+BEARl+Z (oi. 5 )* (RELX1 =0 )* {RELY1<0 ) ) + ( . 5 )* (RELXl =0 
)*{RELY 1>0) 

[28]  BEAR 2‘*( ( " 30PE£y2 *RELX2 ) * ( RELX 2 >0 )x (RELY2iO ))+((( ~30REL 
Y2*RELX2)+Ol)x(RELX2<0))+({ C GORELY2*RELX2 )+02 )x (RELX2 
>0 )*(RELY2<0 ) ) 

[29]  BEAR2+BEAR2+( (oi. S)*(RELX2=0 )*(RELY2<0) )+( . b)x(RELX2=0 
)*(RELY 2>0 ) 

[30]  RVl+({RELXl*2)+RELYl+2)*.5 

[31]  RV2+{(RELX2*2)+RELY2*2)*.b 

[32]  FREQ1+TARZ7]*{ (RV1*TARZGH)* (RELX1*S1XV) * RELY1*S1YV ) * (R 
V^xTARZZ’]  )*  (RELXl  *TARZ3l)*RELYl*TARZ^1 
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6 

f  C  33  3  FRBQ2+TARir\*{{RV2*TARlB])+(RBLX2*S2XV)+RBLY2*S2YV)*(R 

t  V2*TARlB’])  +  (,RELX2*TARlS'})  +  RBLY2*TARlk'} 

\  C34]  L<‘(pTINBN)l  1] 

|  •  [35]  BBAR-(2tL)p(6BARl+L  RNNL  1 + (“2 +NBASNZ 1 t ] ) ) tBEAR2 +L  RUN 

I  L<.U(~2*NEASNi2i]) ) 

[  [36]  FRBQ*-(2tL)p(FREQl+L  RNNL  (  3*V£AS#[  1 ;  ]  ) )  ,FRBQ2+L  RNNL  (31 

'  AfSAS#[2;]> 

[  [37]  SBRX+(2tL)ptSlXPtS2XP 

i  [38]  SBRY*-(2%L)ptSlYPtS2YP 

[  i  [39]  /?SDKX<-(2,L)p,Siyy  .S2W 

f  [40]  flffLyy-(2,L)p,siyv,s2yy 

,  ‘  [41]  DArJAf+l  12p0 

I  [42]  N+ 1 

[43]  s*rixmi!i] 

E  [44]  042,4tf[l|]*l,5.fIJVffJtfCl|2]  ,B£>U?tS;l]  ,Pfffifl[S;l]  .SBMCS;  1  3 

[  ,S5/?y[S;l],J?ff£mS;l].J?ff£mS;l].m[l],W2tl],((S=l)* 

•  S/?Cffl[5])  +  (S=2)*S/?Ctf2[5] 

1  [45]  N+N+ 1 

t  [46]  -(«T>L)/ 0 

[■  [47]  S-riAfffW[iV;l] 

[48]  0Ar*tf«-D/irAMt[l]tf  ,S,TI#£#[0;2]  ,BBA/?[S;W]  .Fif£fl[S;#] 

[S;i»]  ,SB/?y[S;W]  ,BSLVy[SjB]  ,/?Bt^y[S;<n  ,B/1[B]  ,ffy2[B]  ,(( 
S= 1 )xSBCff  l[  5 ] )+ (S=2 )*SRCH2i 5 ] 

;  [49]  -45 

V 
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V00PPFffC2[O]V 

VD0PPVEC2  \M1\M2\M$\R‘%  M  j  DOPPV ;  VAX  \VA1\ NUD ;  f'l  TAM ;  TBMiFAM 
;  FBMi  TDNUBMi  TDNUAM  \  PSI ;  PSID\A  \  \  \ A22 \ A3  3 \ A\ \ A2 \ A3 \ A\2 \ A 
13;<423jBll;B22}Bl  }B2  iB\2\F\REP 

[1]  LOOP* 1 

[2]  Nl+20DATAMt ;4] 

[3]  A#2-*-10i3>ITi4A#C 

[4]  N$+-(*DATAMi  ;  5  ]  )  x  TAR  [  8  ]  +  (Ml  xDATAMi  i&])+M2*DATAMi  ;9] 

[5]  R+ipD AT AM) ill 

[6]  AMP, 3)p0 
[  7  ]  Mi s 1 ]+N 1 

[8]  #[;2>A#2 

[9]  Af[;3]<-A/3 

[10]  P0PPF^-lT/1P[8]x(fflW)+.x(  (B.Dpl) 

[11]  DELTASQ++  /  ,  ( (Af+.  xB0PPF)  +  r-4P[8])*2 

[12]  VAX+DOPPVi 1 S 1 ] 

[13]  IMY«-00PPF[2;1] 

[14]  RUD-DOPPVi 3 ;  1  ] 

[is] 

[16]  W+(D0PPVi  1 ;  1 ]*2 )+00PPY[  2  ;  1  ]*2 

[17]  +(.WiTARiS'}*2)/Cl 

[18]  >I11+-+/W1*2 

[19]  A22++/M2*2 

[20]  i433-«-+/W3*2 

[ 2 1 J  A12++/M\*M2 

[22]  413^-+/Wlx-W3 

[23]  A2Z++/M2x-M3 

[24]  A 1-*-+  /Ml 

[25]  A2*-+  /  M2 

[26]  A3++/-M3 

[27]  Bll^(/113*2)-i411x>J33 

[28]  B2  2-*-(A23*2)-A22*A33 

[29]  Bl-*-r^P[8]x  (/113x/13)-i433x/ll 

[30]  B2*-TAR  [8]*(423Xi43)-433xd2 

[31]  B12-<-(/4  13x^23) -  412x^33 

[32]  F-*-  ( <4  3  3*DELTASQ)  +  (2’i4B[8]*2)x(/43*2)-Bxi433 

[33]  PSI* 3  lpO 

[34]  PS:  PSIilin*-(BllxVAX*2)  +  (B22xVAY*2)  +  (.2xB12xVAXxVAY)  +  (2* 

BlxVAX ) 

[35]  PSI[1 ;1]«-PSI[1;1 ]+( 2xB2xVAY)+F 

[36]  PSI  [2;l]«-(F/lY*2)  +  (ll'/Y*2)-2Mi?[5]*2 

[37]  PSI  l  3;  1  ]■*■(  B1  2  xV  AX  *2  )•*•(  ( B2  2  -B1 1 )  x  VAX  x  V  AY)  +  ( B2  xV  AX )  -  ( Bl  x 
F/1Y)  +  B12xF/iy*2 

[38]  PSID*-  3  3p  0 

[39]  PSIDi  1  ;1>(  2xBllxF/JY)  +  (2xB12xF/4y)+2xfll 

[40]  P5IB[1  ;2]^(  2xB22xF/4y)  +  (2xB12xF^y)+2xfl2 

[41]  PSr0[l;3]«-i43  3 

[42]  PSIDi  2;3]-*-PBIB[3;3]^0 

[43]  P5JB[  2  ;l]«-2xF/1Y 

[44]  PSIDi2',2]*-2xVAY 


*F^™pTT^r™iwi*'i-™»,P!p  r" '  ''HWiflIBPIP 
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rum  p<?rnr3»l1«-(2xBl2xlMX)  +  ((B22-Bll)xVAY)  +  B2 

III]  P!mS;sMt»!2-»u).v«)-(s»m.v«).B1 

[tl] 

[491  -*-R3 

[50]  R2iREP1-RBP-(BPSID)+ ,*PSI 

C513  -*>((  +  /,  [REP1~REP)<.0^)/QS 

[52]  REP- RE PI 

[53]  R3iV AX-RE Pll ;l3 

[55]  F^U33xRBP[3{l])  +  (r4B[8]*2)x(i4  3*2)-i?x^33 

[56]  LOOP-LOOP*  1 

[57]  -PS 

[58]  QS:  VAX-REPi.il  il  J 

[60]  Nt/D*X(  4 13x74X1  +  (A23*V  AX)  +  TAR  [8]x43)*433 

[61]  Cl :STV-NUD ,VAX  tVAX 

V 


[13 

[2] 

[3] 

[H] 


VRMWL[D]V 

VZ-N  RMNL  P'.AAiY 

y*h<  (?Wp  iooooo  >*io  oooi )+.  ooo  oooi 
44^((((l-y)*-if6. 158)-1 >*1*4 . 874) 

AA-(AA- <((y*-l*6.185)-l)*l*4.874))*.323968 

Z«-P[l  ]+44xp[2] 

7 


iSIGMi 03V 

V  SIGMAS Q- SI GM 

[1]  SUVECT  STV 
[23  SIGMASQ-* / (NUAV 

V 


STV 


DATAMl ; 5 ] ) *  2 


x 
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[13 

[23 

[33 

[43 

[53 

[63 

[73 

[83 

[93 

[103 

[113 

[123 

[133 

[143 

[153 

[163 

[173 

[183 

[193 

[203 

[213 


[13 

[23 

[33 

[43 

[53 

[63 

[73 

[83 

[93 


VNUVECTLQlV 

VNUVECT  STV \TB\BB\  TTV ;  PTX V\PTYV 
CVB+DATAMi ; 3  3 <DATAMl ; 1 2  D 
CVA-DATANi  -,3l*DATAMl  ;  1 2 3 
TB+-DATAMI 1 ;  3  3 
BB+DATAMl 1 ;4  3 
TTV<-DATAMi  ;33-TB 

PTXV+ { (DA TAMi ; 2  3  =  1 ) * SRCH 1 [ 3 3 )  +  ( DATAMi ; 2 3  =  2 ) * SRCH2i 3 3 
PTYV+{ { DATAMi }2 3= 1 ) * SRCH 1 [ 4 3 ) + {DATAMI ; 2 3=2 ) * SRCH 2[ 4 3 
RRX+{ (TTVxSTVZ2l-DATAMl { 8  3 ) x CVB)  +  ( ( {DATAMi ; 1 2  3 -TB)*STV 
121- PTXV)  +  (TTV+TB- DATAMI  ;12l)x STVi2l-D AT AMi s 8 ] )xCVA 
RRX+RRX+ {{DATAMi ;2 3 = 1 ) * {RNx 20BB) ) + ( DATAMi ;2 3=2 )x {RN*2o 
BB)  + {SRC  H  lil]+ TB*  SRCH  lit]) -{  SRC  H2ill  +  TBx  SRC  H2ik'}) 
RRY+{{TTV*STVi3l- DATAMi ; 9 3 ) *CVB) + ( ( {DATAMi l 1 2 3- TB) *STV 
[3  3-  PT  XV )  +  {TTV  *TB- DATAMi  ;  12  3  )  *STVi  3  1- DATAMi  ;9  3)*CT>J 
RRY+RRY+{ {DATAMi }2 3 = 1 ) * {RNx lOBS) ) + {DATAMi ;2 3=2 )x {RNxio 
BB)  +  {SRCH  1[  2  3  +  TBx  SRCH  1[  4  3  )  -  ( SRCH 2[  2  ^TBxSRCH 2[  4  3  ) 

RRV+{ {RRX*2)+RRY*2)*.S 
RRX+-RRX-{RRV-RRVil1)xSTVi2l*TARiQ] 

RRY+RRY- { RRV-RRV i 1 3 ) *  STV [ 3  3  *  TAR [ 8  3 
RRN+{ {RRX*2)+RRY*2)*.S 
CO  SV+RRX ♦ RRN 
SINV+RRY*RRN 

CRVS+-TARiS)  +  {  {  {C0SVxDATAMiiBl)  +  SINVxDATAMii9l)xCVB)  +  {{ 
COSVxDATAMi  i8l)  +  SIlfVxDATAMi  ;9  3  )xCVA 
CRVT-TARi9l+{COSVxSTVi2])+SINVxSTVi3] 

NUAV-rSTVi  llxCRVSACRVT 
TTVR4-TTV-RRN*TARi8'} 

V 


VFCtf[03V 

7 FV+FCN  STV ;  EV ;  F 1 ;  F 2 ;  F 3 ;  PRX  TD1  \PRXTD2  \PRYTD1  \PRYTD2 

NUVECT  STV 

EV+NUAV- DATAMi ; 5  3 

F1+-2X  +  /EVxCRVSACRVT 

PRXTD1+{TTVR*RRN*3)x{RRN*2)-RRX*2 

PRYTD l+PRX TD2-{ TTVR* RRN * 3 ) x -RRXx RRY 

PRYTD2+{TTVR*RRN*3)x{RRN*2)-RRY*2 

F2 «-2  x  +  / EV x  { STV [  1  3  ♦  CRV  T* 2  )  *  ( CRV Tx  {D  A TA M [  { 8  3  *  PRX TD\ )  +  DA T 
AM  i  }  9  3  x  PRX  TD  2  )  -  CRV  Sx  CO  SV  +  ( STV  [  2  3  *  PRX  TD 1  )  +  STV  [  3  3  *  PRX  TD  2 
F3+2X+ /EVx{STVil1* CRV T*2)x {CRV Tx {DATAMi i8hxPRYTD\)+D AT 
AMi i$]x  PRYTD2) - CRV Sx SINV+{ STV i 2 3  *  PR YTD 1 )  +  STV [ 3  3  * PRYTD2 
FV+ 3  lp(Fl,F2,F3) 
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VKITERATEm? 

1KITERATE  STV  \II'tCHIS't  LNUB\UNUB\ IVM\FCNB\ STVB\CYC ;  SIGN 
B 

Cl]  II-l 

[2]  NUVECT  STV 

[3]  CHIS4-(pDATAM)in*MEASMlliS]*2 

[4]  LNUB-(DATAMl  ;  5  ]x  ( TAR  C  8  ]- TAR  C  5  ])  tCRVS)  C(  kDATAMZ  ;5]x(r>l/? 
l6]-TARZSH)*CRVS)liyi 

[5]  UNUB-(DATAMZ  ;5  ]x  ( TAR  [  8  ]  +  2MB  [  5  ]  )  *CRV  S)  C  (  ID AT AMI  i*>l*(TAR 
Lfn  +  TARZ5])*CRVS)  Cl]] 

C  6  ]  SET'.  CYC-0 
C  7  ]  SIGHB-SIGM  STV 
C  8 ]  IVN-fti JACOB  STV) 

L9]  FCNB-FCN  STV 
CIO]  STV B- STV 

Cll]  STV A— STV— 3 p ( ( 3  IqSTV ) -IVM+ . xFCNB) 

Cl2]  TEST: R2—(STVZ 2 ]*2 )+STVZ 3 ]*  2 

C 1 3 ]  -(( {STVZll<lNUB)+(STVlll>UNUB)+(R2>TARlS]*2)+( {SIGH  ST 

\ nzSIGH  STVB))<1)/C0N 
Cl4  ]  CYC-CYC+ 1 

CIS]  STV  A— STV— 3  p  (  (  3  lpSfl^B)  -  (  2* -CYC)  x  IVM+  .  xFCNB  ) 

Cl6]  —(.CYCiO  )  / TEST 
C 1 7 ]  STV -STV A- STV B 
C  1 8  ]  CON '.NUVECT  STV 
C  1 9  ]  1 1- 11+ 1 

C20]  SIGMA- SIGN  STV 
C21]  -*((  |  SIGMA- SI GMB)  <10*~5)/0 

C22]  -{II<REPIT) /SET 
V 


V^COBCCHV 

VJA-JACOB  STV  ',CIi  DEL  iJA;FJl  \  FJ2  ;  DEM 
Cl]  Cl- 1 

C  2 ]  DEL— (0,0, 10* -6) 

C  3 ]  JA- 3  OpO 

C  4  ]  FJ1-FCN  STV 

C  5 ]  JI: DEL-~ IfyDEL 
C  6  ]  DEM-DELx( 100  1  l) 

C  7 ]  FJ2-FCN (STV+DEM) 

C  8  ]  JA-JAAFJ2-FJ\)**DEMZCn 
C  9 ]  CI-CI+1 
CIO]  -*-(3  ZCI)/JI 
7 
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VSSSTMZUlV 

VSSSTM  STV ; J ; RR ; K ; AV { RN S H ; STV W ; FM 

[1]  «/«-l 

[2]  SSTM+ 32  6pO 

[3]  SRNiRN-10*J 

[4]  KITE RATE  STV 

[5]  SSTMZ Jlh-RN .STV A  .SIGMA. J 

[6]  t 1—J  + 1 

[7]  -(JS12  )/SRN 

C 8  3  AV- 12p0 

[9]  dK[  1  ]—SSTMl  1  ;5  ]< SSTMZ  2  {5  ] 

[10]  AVI.  12  1-SSTMZ  12  ; 5  SSTMZ  11  ;5  ] 

[11]  K-2 

[12]  VECiAVtm-{SSTMZKlS']<SSTMtK-l  ;5  ]) 

[13]  AVlK]-AVZn*SSTMZKiSl<SSTMiK  +  l  ;5] 

[14]  K-«-K+l 

[15]  -*-U<12)/VEC 

[16]  -( (+/ AV)*1)/ SUP 

[17]  SSTMZ  1 3  ;  ]«-6  p  0 

[18]  -*-0 

[19]  SUP: RN  S-AV / 12  +  .SSTMZ ;6 ] 

[20]  N- 1 

[21]  3rr:ffiV«-SSrM[J?WS[An;l]-5 

[22]  S2TV-SSrM[iWS[ff]  }2  3  4] 

[23]  KITERATE  STVU 

[24]  SSTMZ  12  +  {  3*N) -2  ll-RN  tSTV  A, SIGN  A,  12  + -2 

[25]  RN-SSTMZRNSZtni  l] 

[26]  KITERATE  STVW 

[27]  SSTMZl2+(3*N)-l  ;]-RN , STV A .SIGM A . 12+ (  3xff)-l 

[28]  RN-SSTMZRNSZN1 Si ]  +  5 

[29]  KITERATE  STVW 

[30]  SSTMZ  12  +  3  xff;  ]*-/?fl  ,S2T>I  .SIGMA  ,  12+3*ff 

[31]  AKW+1 

[32]  -(N&pRNS) /STT 

[33]  /Mfc(12  +  3xpfftfS)  +  ,SSm  ;5])[1] 

[34]  RN-SSTMZRi 1 ]-2 . 5 

[35]  STVW'-SSTMZR’,2  3  4] 

[36]  KITERATE  STVW 

[37]  FM- 3  5p  0 

[38]  F.V[1  ’.1-RN  .STVA.SIGMA 

[39]  RN-SSTMZR 1 1 ] 

[40]  KITERATE  STVW 

[41]  FMZ2-.1-RN  .STVA.SIGMA 

[42]  RN-SSTMZRin+2  .  5 

[43]  KITERATE  STVW 

[44]  FMZ3l1-RN .STVA .SIGMA 

[45]  MIN—FMZ  i  h.  FMZ  ;5])[1]{] 

V 
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VRRNG [0]7 

7  RRNG  STV ;  NUA  ;  TUA  ;  BB\RVV ;  EVV ;  RIM ;  RIMD  }  RCN  { DRCN ;  RNA ;  DN-, 
DPDiLBiFViUB 

C  1  ]  NUA-iQtDATAML  ;2  3  =  J  )x  (DATAMl  } 3  USRCHll  5 ]  )  ) /<*DATAMl  ;5  ] 

[2]  TUA+((*(.DATANl  {2  ]  =  1  )x  (DATAMl  i  3  IsSRCH  lC  5  ]  )  )/*DATANl  ;3]) 
-DATAMlUSl 

[3]  BB+DATAM [l;4] 

[4]  RVV-STV [  2  3  3 - 5/PC* » 1 C  3  4] 

[  5  ]  EVV+NUA-STVZH 

[6]  DPD-{  (  20BB )  *RVVl  1  ]  )♦  (  loflfl)  x/?yy [  2  ] 

[7]  (( </?mi3*2)  +  fll't'[2]*2)*.5)x(pfftfy|)-i)  +  0PD 

[8]  □■»-F,^(r>4i?C8]*Sr^Cl])x  +  /^7 

[9]  0<-UB+(pNUA)x-DPD 

[10]  +(UB&FV)*(FV&UB)  =  1)/C0NT 

[11]  -*{FV<LB)  / LBL 

[12]  +{UB<FV) /UBL 

[13]  LBLiRN+10 

[14]  UBLiRN* 120 

[15]  *  MO  INITIAL  RANGE  SOLUTION  POSSIBLE ,  RN  =  *  ;/W;  ’  ASSU 
MED  * 

[16]  -*0 

[17]  CONT:RIM+(2,(pNUA))pO 

[18]  RN+  50 

[19]  CALC : RIMD+RN+ TUA x DPD 

[20]  RIMZ  1{X  ( RNx20BB)+TUAxRVVZll)*RIMD 

[21]  RIMZ2il4-{(RNxlOBB)+TUAxRVVZ2l)*RIMD 

[22]  /?Cff^((Tiifl[8]*Smi])x  +  /£m  +  (/?mi]x+/fl.rv[l;])  +  /?m2]x 
+/RIMZ2i1 

[23]  DRCN--(((RVVZll*2)+RVVZ2'}*2)-(DPD*2)  ) x+ / TUA* RIMD+2 

[24]  RNA+-RN-  RCN*DRCN 

[25]  DN*-\RNA-RN 

[26]  RN+RNA 

[27]  +(DN>10*~* )/CALC 

[28]  ' NEWTON- RAPH SON  SOLN  FOR  RN  =  • iRN 
7 


7flffSff£[[]37 
7 RNSEL  STViC 

[1]  RNMT+23  2p0 

[2]  O 1 

[3]  MS:ff^l0  +  5x(C-l) 

[4]  RNMTiCi >JW ,SIGM  STV 

[5]  C«-C+l 

[6]  -►(CS23  )/RNS 

[7]  RN^RNMTZ(b,RNMTZi2l)Llhll 
7 
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vsffimcnv 

7  SHIFT 

Cl]  * SET  XTV  tYTV  tXSV  tYSV' 

[2]  iy«-,o 

[3]  X2V«-ini] 

[4]  YTV+IVW 

[5]  X57«-IV[3] 

[6]  7SF^I^[4] 

[7]  ' SET  H  =  NO.  RBPIT. • 

c  8  ]  #«-,  n 

[9]  • SET  D  *  IHTBR  PULSE  TIMS' 

[10]  D+t □ 

C 1 1  ]  R+ 0 

[12]  C’ff5:X«-?150 

[13]  X^?150 

[14]  '(/HEN  X  -  ’;X;'  AND  Y  =  ‘{X 

[15]  DPXACr 

[16]  0APPP0X 

[17]  'DBXACT  =  •  ;DA1 

[18]  » DAPPROX  =  » ;0AA;‘  PERCENT  ERROR  =  •  ;  1  OQx(.DAA-DAl  )*DA1 

[19]  'DAPPROX2  ~  '  iDAA  2;’  PERCENT  ERROR  =  • ;100x(DAA2-DAl)* 
DAI 

[20]  fr*-X+l 

[21]  K&N)  /  CHS 
7 


7DEXA02'[D]7 

VDEXACT 

[1]  0A«-D 

[2]  DAl^D+(*C)x(  ((  (X-(DAxXSy)-DxXrt')*2)  +  (r-(DAxy57)-Dxyri') 
*2)*. 5)-( (X*2)+y*2)*.  5 

[3]  ■+((  lDA-DAl)<10*~10)/0 

[4]  DA+DA 1 

[5]  -2 
7 


7DAPPP0X[D]7 
ID  APPROX 

[1]  fl«-((X*2)+y*2)*.5 

[2]  DAA-*-Dx(C+(X2VxX»P)+yrKxy*P)*C,+  (X5FxX*P)+yS^xy*D 

[3]  0AA2«-Dxl  +  (  ♦C,)x(X2,7xX*P)  +  (yrFxy*P)-(  (XSV*X*R)+YSV*Y*R) 
7 
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NSHIFTYlUlV 
7 SHIFTY 

Cl]  'SET  XTV  tYTV  tXSV ,YSV  t  X.Y' 

C  2  ]  IV«-  ,□ 

[3]  X2V«-IVCl] 

[4]  YTV+IV12] 

C  5  ]  X5^IK[3] 

[6]  yS^I^C4] 

[7]  X«-IVC5] 

[a]  x^-me] 

[9]  *  SET  N  =  NO .  REPIT. * 

CIO]  #«-,  D 

Cll ]  1  SET  D  =  INTER  PULSE  TINE' 

C 12  ]  0«-,O 

C 13  ]  X«-0 

C 14 ]  CHSiX+X+.l*XTV-XSV 

c  1 5  ]  r+r+.  lxrry-xsy 

C 1 6 ]  'WHEN  X  =  'iXi'  AND  Y  =  '  ;Y 
C 1 7  ]  DEXACT 
C 1 8 ]  D APPROX 
C  1 9 ]  'DEXACT  =  »  -,DA1 

C  20 ]  ' DAPPROX  =  'iDAA;'  PERCENT  ERROR  =  * ; 1 00  * (DAA-DAl ) *DA 

C 21 ]  'DAPPROX 2  =  *  \DAA2\ *  PERCENT  ERROR  =  » { 1 OOx (DAA2-DA1 ) 
DAI 

C  22  ]  fr-K+1 
C  23 ]  -*• (KiN)tCHS 
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APPENDIX  B 

THE  UNIQUENESS  PROBLEM  FOR  PERFECT  INFORMATION 

If  there  were  no  errors  in  the  measurements  of  bearing  or 
frequency,  the  fundamental  system  [B-2],  which  derives  from  several 
measurements  of  [B-1],  would  possess  a  unicue  solution  for  all 
parameters,  provided  at  least  three  different  bearings  were  included. 

fa  (C  +*r  .  *VT)  *  f  (C  +*r  .  %)  IB-1] 


(C+?1 

<C  +  TH  ■  \>  f  -  \  cos  BN  *1 


sin 


f 

*N 


sin  bn  Yt 


[B-2] 


To  prove  the  assertion  of  a  uniaue  solution  for  the  ranee 
parameter,  it  suffices  to  establish  that  the  system  [B-2]  has  a  uniaue 
solution  for  the  three  parameters  f,  XT,  Y,p.  Indeed,  Geometrical 
considerations  make  it  clear  that  any  two  distinct  bearinas,  and  B^ , 
obtained  at  times  TT^  and  TT^  (without  error),  together  with  the  known 
values  XT>  YT  will,  to  all  practical  purposes,  uniauely  determine  ranae, 
the  very  small  Indeterminacy  introduced  by  the  change  of  transmission 
time  of  the  two  pulses  being  nealiaible. 


To  prove  that  b-2  possesses  a  uniaue  solution,  we  only  have  to 
consider  N=3  with  B^,  Eg  and  B^  all  different.  Under  these 
circumstances  there  will  be  a  uniaue  solution  so  Iona  as  the  determinant 
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D  « 


C  +  r. 


C  +  r. 


C  +  r. 


-  f  cos  B1 

-  fg  cos  Bg 

-  f  cos  Bg 


-  f  sin  B1 

-  f2  sin  Bg 

-  f„  sin  B. 

o  3 


[B-3] 


does  not  vanish. 

Eauation  [B-1 ]  makes  it  possible  to  replace  the  elements  of  the 
first  column  of  [B-3 3  by  expressions  of  the  form: 


(C  ^ 


Then  the  determinant  can  be  cleared  of  the  f^  by  multiplying  each  one  of 
its  rows  by  the  appropriate  non-zero  factor  f/f^.  The  following 
simplified  determinant  D'  thus  results.  D'  =  0  if,  and  only  if,  D  =  0. 


D* 


C  +  cos  B1  X^,  +  sin  B1  YT 
«  • 

C  +  cos  Bg  +  sin  Bg  YT 
•  • 

C  +  cos  B^  Xjp  +  sin  Bg  Y,^ 


-  f  cos  B^ 

-  f  cos  Bg 

-  f  cos  Bg 


-  f  sin  B1 

-  f  sin  Bg 

-  f  sin  B^ 


[B-L  ] 


The  determinant  D'  can  vanish  only  if  there  are  constants  "a”  and 
"b",  so  that  the  third  row  is  the  linear  sum  of  "a"  times  the  first  row 
and  "b"  times  the  second.  Then,  the  three  equations  in  [B-53  would  be 


satisfied : 
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Bin  ■  a  sin  +  b  sin  Bg 
cos  B^  ■  a  cos  B1  +  b  cos  Bg 

C  +  cos  B3  Xt  +  sin  B3  Yt  -  (a  +  b)  ^  [B-5] 

C  +  (a  cos  B1  +  b  cos  Bg)  X^, 

+  (a  sin  B1  +  b  sin  Bg)  YT 

However,  substituting  the  first  two  equations  into  the  last  yields 

C  s  (a+b)C  or  a+b  =  1  [B-6] 

Simultaneously,  the  first  two  equations  of  [B-5]  state  that  the  unit 
vector "r*3  is  the  linear  combination  ar^  +  br^.  However,  vectors  of  the 
form  a?1  +  br2>  with  a  +  b  =  1,  stretch  from  the  origin  in 
two-dimensional  space  to  points  on  the  line  joining  the  termini  of  the 
two  unit  vectors  r1  and  r  •  However,  it  can  be  seen  geometrically  that 
no  other  vector  to  that  line  has  length  one,  except  r1  and  rg 
themselves.  Thus  since  no  two  of  r^,  r?,  r3  can  coincide,  the 

determinant  D  cannot  vanish,  and  the  expressions  in  [B-3]  yield  unioue 
solutions  for  f,  XT,  YT,  as  assert'1 
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